Density waves and Cooper pairing on the honeycomb lattice 
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Motivated by the surge in research activities on graphene, we investigate instabilities of electrons 
on the honeycomb lattice, interacting by onsite and nearest-neighbor terms, using a renormalization 
group scheme. Near half band-filling, critical minimal interaction strengths are required for instabil- 
ities toward antiferromagnetic or charge-density wave order. Away from half filling, /-wave triplet 
pairing and d + id singlet pairing instabilities are found to emerge out of density-wave regimes. 



Recently it has become possible to fabricate 
graphene PJ. Already the structure of graphene poses 
some interesting questions to experiment and theory 
Concerning the electronic properties, the massless-Dirac 
spectrum near the two Fermi points of the undoped 
system gives rise to novel quantum Hall effects and 
many more complex physical puzzles [2]. Unconven- 
tional mesoscopic properties of graphene as specular 
Andreev reflection [3[, analogues to relativistic quantum 
mechanics Q or spin quantum bits with graphene [BJ are 
likewise discussed intensively. 

Many-body effects are expected to play a role in 
graphene, although no drastic anomalies have been re- 
ported yetQ. Superconductivity in graphite compounds 
is usually attributed to extrinsic causes, e.g. to inter- 
layer states for intercalated compounds[7j or, more in- 
triguing and less clear-cut, to disorder with sulfur atoms 
in graphite-sulfur composites Q. In theory, the elec- 
tronic selfenergy for the undoped case was found to be 
of marginal Fermi liquid type[9| . In the doped case, 
a normal Fermi liquid develops To[ . The Dirac spec- 
trum should further be reflected in the dependence of the 
phonon renormalization on the electron density 11| . Re- 
cently it has been proposed that a strong second-nearest- 
neighbor repulsion may lead to a nontrivial insulating 
state with a quantized spin Hall conductance (12J. There 
have been various theoretical studies of possible ordered 
phases on the honeycomb lattice ID, 14 , TBI, III 17 1 driven 
by interactions, but to date no comprehensive picture 
for a wider range of interaction parameters and band 
fillings is available. Here we investigate the intrinsic 
interaction effects of electrons on a honeycomb lattice 
for short-ranged interactions, using a perturbative func- 
tional renormalization group (fRG) technique. The fRG 
is known to give an unbiased picture of the leading or- 
dering tendencies of interacting fermions at low temper- 
atures. Recently it has been applied to two-dimensional 
(2D) Hubbard models on square[18|, |l{J, [2fJ, |2l| and tri- 
angular lattices 22 , 23| . 



Our model is a 2D honeycomb lattice with nearest- 
neighbor hopping amplitude t. The interaction terms 
contain onsite and nearest-neighbor repulsions U and V, 
and a spin-spin interaction J. The restriction to short- 
range terms is partly due to the difficulty to treat a long- 



range Coulomb part directly in the fRG approach. How- 
ever, most experimental graphene systems are doped to 
some degree. Hence at least the effective interaction is 
screened, and our starting point may not be unrealistic. 
Furthermore, for the undoped case, the long-range part 
was shown [15j to be marginally irrelevant in 1/N, hence 
our results for short-range interactions may even be use- 
ful in the undoped case. The Hamiltonian reads 
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denotes all pairs of neighbored sites and $$ = 
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The fRG scheme used here is an approximation to an 
exact flow equation for the one-particle irreducible ver- 
tex functions of a many-fermion system when a param- 
eter in the quadratic part of the action is varied [24j. I n 
the temperature-flow scheme employed here, the temper- 
ature T is used as flow parameter. The fRG flow is gen- 
erated by lowering T from an initial value To where the 
interaction effects are negligible. In the approximation 
we use 2l| , the change of the interaction vertex Vt with 



is given by the T-derivative of one-loop particle-hole di- 
agrams, including vertex corrections and screening, and 
particle-particle diagrams, of second order in the vertices 
Vt- Higher loop contributions are generated by the inte- 
gration of the flow. Like in many previous works using 
this method, the self-energy feedback on the flow of Vt 
is neglected. It may become important when the inter- 
actions get large. Hence the flow is stopped when the 
interaction strength exceeds twice the bandwidth. The 
method is controlled in the limit of interaction strength 
going to zero. For the interesting case of moderate inter- 
actions, it should be viewed as a step beyond the mean- 
field level that captures the evolution and competition of 
various correlations in an unbiased way. 

The interaction vertex can be expressed by a coupling 
function Vy(fci, &2, ^3). It depends on the generalized 
wavevectors of two incoming particles (fci and fe) and one 
outgoing (£3) particle with wavevector, Matsubara fre- 
quency and spin projection fcj = (fcj, u>i, Si). In the search 
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for instabilities toward symmetry breaking, the frequency 
dependence is neglected and the u>i are set to zero. The 
/c-dependence of remaining function Vr{ki, k2,k$) is dis- 
cretized in the so-called TV-patch scheme, introduced 
in this context by Zanchi and Schulz[l8j]. For stan- 
dard many-fermion systems with a FS, this amounts to 
keeping Vr(ki, U2, £3) constant within patches labeled by 
k{ = 1, ...N perpendicular to the FS. This defines an 
iV 3 -component coupling function Vr(klk2, k^), which is 
computed for k(ki) = hi, i = 1, . . . 3, on the FS. As for 
the honeycomb lattice there is no FS at half band filling, 
we generalize the patching scheme to two or three rings 
of 18 or 24 patches around the Dirac points (see Fig. []}. 
In addition we have a band index for the incoming and 
outgoing particles. 

The RG flow is started at an initial temperature To. 
The initial interaction is given by the bare interaction 
with the onsite repulsion U and the nearest-neighbor re- 
pulsion V. Specifically, we search for flows to strong cou- 
pling, where for a certain low temperature T c one or sev- 
eral components of Vr(fci,fc2,&3) become large. At this 
point the approximations break down, and the flow has to 
be stopped. Information on low-temperature state is ob- 
tained by analyzing which coupling functions grow most 
strongly and from the flow of susceptibilities. In partic- 
ular for 2D systems, this instability does not guarantee 
true long-range order. Rather, it should be interpreted as 
a breakdown of the (semi-)metallic state and as indicator 
for the leading correlations at low T. 

We begin with the semimetal for zero doping and chem- 
ical potential /i = 0. The first finding is perturbative 
stability. Starting the RG flow at high T and small U, 

V or J, we can follow the flow down to lowest T without 
a divergence. This is quite different from many-fermion 
systems with a finite density of states at the Fermi en- 
ergy, where the flow practically always leads to some kind 
of instability. In our case, the absence of a flow to strong 
coupling indicates the absence of long-range order due to 
electronic interactions even at lowest T. 

Next we increase the bare interactions. Above a crit- 
ical value U c ~ 3.8t for the onsite repulsion U with 

V = J = 0, the interactions flow to strong coupling. The 
static antiferromagnetic (AF) spin susceptibility grows 
most strongly toward the critical temperature scale, in- 
dicating a tendency toward AF spin-density wave (SDW) 
formation with opposite orientation of the ordered spin 
moment on A and B sublattices. For small U, J = 0, and 
increasing nearest-neighbor repulsion V we again find a 
flow to strong coupling for V > V c ~ 1.2t, now with lead- 
ing charge-density wave (CDW) correlations for different 
charge densities on the two sublattices. This compares 
favorably with a previous l/N- analysis [l5| finding the 
same instabilities beyond critical values U c and V c . There 
is also good agreement with the early Quantum Monte 
Carlo work by Sorella[l3j who found a transition to an 
AF Mott-state at U ~ 4i. In Fig. Q] we show the depen- 



dence of the critical temperature T c for the flow to strong 
coupling on the interaction parameters. The critical U 
for boundary of the SDW instability is more or less unaf- 
fected by an increasing V, while V c for the CDW regime 
is shifts to larger V with a roughly linear dependence on 
U. When the two lines meet, there is a continuous change 
in the flow from leading SDW to leading CDW correla- 
tions (or vice versa). This is consistent with a first order 
transition, if order is possible at all. The competition for 
the low energy spectral weight not included in this study 
as the selfenergy is neglected could however reduce the 
ordered moments in the transition region. 

In Fig. [2] we display the flow of various susceptibilities 
for SDW and CDW instabilities, and Vt very close to T c . 
The actual calculation takes place in the fermionic basis 
which diagonalizes the hopping term. The resulting in- 
teractions are transformed back into the sublattice basis 
with operators c~ on sublattice b — A,B for incoming 

k,s,b 

and outgoing particles. In Fig. [2]c), d) we show the real 
part of the effective Vr(fci, &2, fe) for a small number of 
N = 12 fc-space points near the Dirac points (with points 
1,2 and K, 3 and 4 at K' , continuing clockwise around 
the Brillouin zone hexagon) close to T c of the SDW in- 
stability. Patch indices 1 to 12 belong to particles on 
sublattice A, and 13 to 24 to sublattice B. As function of 
the incoming k\, ki with outgoing k% fixed, Vr(ki, fe, k%) 
shows either vertical or horizontal features with strongly 
attractive or repulsive values. The vertical features have 
k% = &3 (i.e. same wavevector, same sublattice for parti- 
cles 2 and 3) or &2 = &3 ± 12 (same wavevector, but dif- 
ferent sublattices for particles 2 and 3). We can compare 
Vr(ki, &2, ^3) with an infinite-range interaction which 
gives a SDW as groundstate. On a lattice with TV sites, 
we can define spin-spin interactions iV -1 J2qJ$ b 
with S~=h y~*r cr , 1,. For the infinite-range 

q 2 ss k+q,s,b K,8',b to 

SDW interaction, only the q = components in J^b,y 
are nonzero. We should have JgL < 0, i.e. ferromag- 
netic (FM) on the same sublattice, while J-4 Q > 0, i.e. 
AF, for different sublattices. Comparing this with the 
effective interaction from the fRG, 

— V T b2b3 (£ k', ^+9)4+g-, S ,6 3 4'-?-, S '^4 C fe'^'^2 C fe, S ,f.l ' 

we get V£ b ' b ' '(k,k',k' - q) = -jb b ' and X^ b ' h (£, k', k + 

q) = -Jp b ' /2. In the fRG data in Fig. Q] c) and d), 
only the q = interactions grow strongly, and the signs 
depending on the sublattice follow exactly that of the re- 
duced spin-spin-interaction with FM intra-sublattice and 
AF inter-subband processes. The CDW instability can be 
read off from Vt in a similar way. 

Next we turn to the doped system. Moving fi away 
from zero, we obtain two FSs around the Dirac points 
and the q = nesting between the two bands is re- 
duced. This cuts off the CDW and SDW instabilities 



3 





k /a 



2 2 
U/t v/t 



FIG. 1: (color online). Left: Brillouin zone and 18x3 points 
used for discretizing the wave-vector dependent interaction. 
The solid lines are at constant band energy. The lattice con- 
stant (minimal distance between two A-sublattice sites) is set 
to unity. Right: Critical temperature T c for the flow to strong 
coupling vs. interaction parameters U and V at half filling 
(i — 0. In the region with T c = 0, the semimetal is stable. For 
small U and V > 1.2t, the flow is toward a CDW instability, 
for small V and U > 3.8t toward a SDW instability. 
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FIG. 2: (color online). Upper plots: Flow of susceptibilities, 
solid line CDW, dotted line SDW, for a) U=4t, V=0, b) (7=0, 
V=1.5t. Lower plots: Effective interactions Vr(ki,k2,k3) 
very close to the SDW instability at T c w 0.05i for U=4t, 
V=0. The colorbar indicates the values of the couplings. The 
incoming wavevectors fci and ki are on 12 points on the inner 
rings near the Dirac points, for k\i% =\ to 12 on sublattice A, 
and fcj/2 =13 to 24 for sublattice B. The 1st outgoing parti- 
cle kz is at point 1 and sublattice A. In c), the 2nd outgoing 
particle is on sublattice A, in d) on sublattice B. 



at low T. The behavior of the CDW susceptibility and 
T c vs. /Lt is shown in Fig. [3] a), b). The SDW insta- 
bility for dominant U (or J) behaves analogously. Be- 
yond a critical doping, the CDW susceptibility remains 
finite for T — ► 0. If we continue the flow down to lower 
temperatures T < 10~ 3 £, we observe a strong growth 
in the Cooper pairing processes with zero total incom- 
ing wavevector. This is clearly visible in the effective 
interactions near T c shown in Fig. [3J c) for /i = 0.75i. 



FIG. 3: (color online), a) Flow of pairing susceptibilities in /- 
wave (solid lines) and CDW channel (dashed lines) for chem- 
ical potentials fj, = O.lt (thick lines) to fi = 0.75t (thin lines), 
b) Critical temperatures for the flow to strong coupling vs. /i. 
Crosses: U=1.2t, V=2At with a CDW instability for /j, < 0.7t 
and a triplet Cooper instability for /i > 0.7t. Circles: U=1.2t, 
J=2At with a SDW instability for /i < 0.65t and a singlet d- 
wave Cooper instability for /i > 0.65t. c) Effective interaction 
at low scales for U = 1.2t, V = 2.4t, n = 0.75t, and outgoing 
wavevectors fixed on sublattice A, k^ — 2 for a 18 x 3 dis- 
cretization. Points fci and ki are on the middle of the three 
rings nearest to the FS and on sublattice A for index 1 to 
18 and on B for 19 to 36. The colorbar indicates the values 
of the couplings, d) Crosses: Pair scattering Vr(fei, fe) 
for [7=1. 2t, V=2At, /i = 0.75t with incoming and outgoing 
particles on sublattice A and ki + fc(fei) = vs. fe around 
the Brillouin zone hexagon. The dotted line is the ansatz 
— Vfd*f(ki)df(k2). Circles: Same data for the d-wave instabil- 
ity at U = 1.2t, J = 2.4i, n — 0.75t. Here the dashed line is 
-V d [d* x2 _ y2 (k 1 )d x 2_ y 2(k 2 ) + dZ y (ki)d X y(k 2 )]. 



Here, processes with zero total incoming wavevector (di- 
agonal features) are enhanced strongly. The pair scatter- 
ing Vr(k, —k — > k', —k') is odd with respect to reversal 
of the outgoing (or incoming) wavevectors, correspond- 
ing to triplet pairing. From Fig. [3] c) it can inferred 
that the pair partners (fc, —k) are on the same sublat- 
tice. The pair scattering between wavevectors near the 
same Dirac point is attractive, and from one to the other 
Dirac point it is repulsive. As shown in Fig[3jd), the zero- 
total-momentum part of the effective interaction follows 
closely the form —Vfdf(k)df(k') with the /-wave form 
factor df(k) = sin(fca,) — 2 sin(fc x /2) cos(y/3k y /2). The 
pairing has the same sign on the two sublattices. The cor- 
responding mcanficld picture gives a nodeless state with 
gap amplitudes of opposite sign on the two Fermi circles. 
In real space, the pairing with df(k) takes place between 
a given site and its 6 next-nearest neighbors, with sign 
change upon a 7t/3- rotation around the site. Thus the 
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pairing is not directly mediated by the nearest-neighbor 
repulsion V, but it is due to a next-nearest neighbor at- 
traction that is generated by second and higher orders of 
V summed up in the fRG flow. Consistently, if we intro- 
duce a next-nearest neighbor attraction V% < 0, the pair- 
ing instability is enhanced. Note that the intra-sublatticc 
pairing found here is quite different from the nearest- 
neighbor inter-sublattice pairing found in a meanficld 
study in Ref. for attractive V. In full agreement 
with the analysis of the interactions, the flow of the /- 
wave susceptibility with form factor df(k) (see Fig.[3]a)) 
shows a strong upturn once the CDW susceptibility is 
cut off by sufficient doping. Other pairing susceptibili- 
ties grow much more weakly. 

The SDW regime driven by the onsite-C/ does not pro- 
duce any measurable scale for a superconducting insta- 
bility within our numerical precision. Upon doping, the 
SDW regime just gives way to a stable Fermi liquid. 
The SDW regime can also be generated with a antifer- 
romagnetic Heisenberg interaction J on nearest neigh- 
bors. The semimetal becomes unstable with respect to 
the SDW for J > 2t. Doping into this state (see Fig. [3] 
b)), the SDW scale breaks down. Now we find an insta- 
bility in the singlet pairing channel with d-wave symme- 
try, that can also be found from a meanfield decoupling 
of the J-term[l7|. The energetically best state is the 
complex linear combination of the two degenerate ba- 
sis function d xy and d x 2_ y 2, leading to a fully gapped 
time-reversal-symmetry-breaking state 17[ . In Fig. [3] d) 
we show the effective pair scattering. It follows closely 
the ansatz -V d [dl^_ y2 {k)d x 2_ y 2{k')+d% y {k)d xy (k')] with 
the form factors corresponding to next-nearest neighbor 
pairing of d x 2_ y 2- and d^-type, d x 2_ y 2(k) — e~ lky /^ — 
e ik y /(2V3) cos k x /2 and d xy {k) = i e lk y/^) sin/^/2. 

Regarding possible realizations of Cooper pairing, we 
note that a sizable T c for the triplet pairing requires dom- 
inant nearest-neighbor repulsion V > U large enough to 
be at least close to a CDW ordered state for zero dop- 
ing. V > U could be realized due to Holstein phonons, 
reducing the effective U 25]. A CDW state (or an an- 
tiferromagnetic insulator indicative of strong Heisenberg 
exchange J which could lead to robust d + id pairing), 
however, does seem to be realized in graphene. Yet, 
for graphene on a substrate, there have been theoretical 
ideas [261] that out-of-plane vibrations of the carbon atoms 
with opposite amplitude on the two sublattices could 
trigger a CDW instability, at least in a magnetic field. 
If the system is near such an instability, the phononic 
fluctuations would add to the intrinsic fluctuations, and 
move the system closer to the parameter range for triplet 
pairing. Recent photoemission work[27| for graphene on 
SiC substrates revealed a gap-like feature near the Dirac 
points, interpreted as AB-sublattice symmetry-breaking 
by the substrate. If this gap is in fact a cooperative 
effect of substrate and electronic interactions, it would 



seem promising to dope the system out of gapped phase 
and to search for superconducting correlations. 

In conclusion, we have analyzed instabilities of inter- 
acting electrons on the honeycomb lattice using a fRG 
method. The undoped state becomes unstable with re- 
spect to spin-density wave and charge-density wave in- 
stabilities, if onsite or nearest-neighbor repulsions exceed 
critical values. Upon sufficient doping, the CDW in- 
stability gives way to a triplet-pairing instability with 
intra-sublattice Cooper pairing of next-nearest neigh- 
bors. Doping of a SDW regime with J > leads to a 
singlet-pairing instability in the d-wave channel. 

CH thanks A. Black-Schaffer, G.H. Gweon and M. 
Kinza for discussions and BaCaTec for financial support. 
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